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Abstract 

We present theoretical predictions of granular flow in a conical hopper based 
on a continuum theory employing a recently-developed constitutive model 
with microstructure evolution by Sun and Sundaresan pp. The model is 
developed for strain rate-independent granular flows. The closures for the 
pressure and the macroscopic friction coefficient are linked to microstructure 
through evolution equations for the coordination number and fabric. The 
material constants in the model are functions of particle-level properties. A 
salient prediction is the variable stress ratio along the flow direction, in con- 
trast to the constant ratio employed in some widely-used plasticity theories, 
but supported by results obtained from discrete element simulations. The 
model permits direct interrogation of the influence of particle-particle friction 
as well as normal-stress differences on the stress distribution and discharge 
rate. Increasing particle friction leads to higher stress ratios, but lower nor- 
mal stress and flow rates, while considering normal-stress differences results 
in the opposite effects. 
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1. Introduction 

The flow of granular materials in bins and hoppers is of great practical 
interest in connection with the handling and transportation of bulk solids 
such as sand, coal, ore and grain. It has been extensively studied for sev- 
eral decades and generated a significant number of publications with earlier 
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results summarized in monographs [21 [3], and review articles [HE]. With ad- 
vances in experimental techniques and computational power, recent results 
using particle image velocimetry (PIV) [HI [TJ El [91 [10] and discrete element 
method (DEM) pH [121 H31 DEI [161 E] have shed light on flow kinematics 
and on particle-property effects on the flow behavior. The understanding 
also drives the development of continuum theories [TBI EE 120], which are 
necessary for predicting industry-scale hopper flow. 

In spite of the extensive research, a continuum model capable of accu- 
rately predicting realistic hopper-flow behavior is still under development. 
For example, the challenges faced in the constitutive modelling of hopper 
flow were described in [21]; deficiencies of a Mohr-Coulomb plasticity theory 
were summarized in [18] . The deficiencies of current models are also reflected 
in the fact that widely different predictions of well-defined silo flow were made 
using different models [22] . Precisely due to these difficulties, the hopper-flow 
problem also serves as a stringent benchmark case to test newly-developed 
theories for granular materials, such as the stochastic flow-rule [IB] , non-local 
hypoplasticity [19] and elasto-plastic models [20] . 

In this paper, we employ a new constitutive model with microstructure 
evolution pQ to predict flow in a conical hopper. We utilize the radial solu- 
tion, first presented by Jenike |23j , to reduce the number of spatial variables 
to one by the inherent symmetry. Such a solution, albeit being rather sim- 
plified, helps to elucidate the essential capabilities of the constitutive model 
and prevents obscuring the features by complexities in boundary conditions 
and numerical methods. We will demonstrate that the model is capable of 
predicting a varying stress ratio along the flow path, which is realistic but 
elusive in previous theories. We will show how the particle properties and 
microstructure are linked to the continuum variables, which enables study- 
ing the role of such microscopic parameters at a much lower computational 
cost. We will also compare the predictions to DEM simulation results and 
the so-called hour-glass theory (HGT) [2H [23 [3J that is also a radial-flow 
solution, but employing the Mohr-Coulomb plasticity theory. 

The remaining part is organized as follows; the hopper-flow model is 
described in section [2] detailing the governing equations (2.1) and the con- 
stitutive model ( 2.2[ ). The results from solving the model are presented and 
analyzed in section [3j leading to the conclusions and future work in section [4j 
Regarding notation, we employ lightface italics for scalars, the boldface reg- 
ular font for vectors and the San Serif font for the second-order tensors. 
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(a) (b) 

Figure 1: The schematic diagram of the geometry and the spherical coordinate systems 
for a conical hopper showing (a) a three-dimensional view and (b) a r-9 cross-section view 
with arrows pointing to the virtual apex indicating the radial flow. 



2. Hopper flow modeling 

2.1. Flow conditions and governing equations 

We consider an incompressible cohesionless granular material of mass 
density p and volume fraction discharging from a narrow-angled conical 
hopper of half angle a, whose geometry is sketched in Fig. [1} We use a 
spherical coordinate system ((r, 9, ip) see Fig. |I{a)) with origin at the virtual 
apex of the hopper. The walls are assumed to be rigid and smooth. Since the 
angle a is small, gravity is assumed to be in the radial direction. By axial 
symmetry, the shear stress a r g is reasonably assumed to be zero. We further 
assume that the velocity is radial and independent of 9 and (p coordinates. 

With these conditions, the continuity equation for this incompressible 
granular flow is reduced to 

V 



V r = — - 



(1) 



where V is a positive flow rate independent of r. This flow rate is related to 
the discharge mass flow rate W, as W = 2np(j)(l — cosa)V. We then obtain 
a dimensionless momentum conservation equation utilizing Eq. [T] 
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where r is scaled by the radius of a free-fall arch r (see its location in 
Fig. [T | (b ) ) , i.e. f = r/rn\ and the normal stresses a = a/(pgr ) and V = 
V I \/r\g. Solution of Eq. ^ with appropriate boundary conditions predicts the 
radial distributions of velocity and stress as well as the flow rate. However, 
this equation is not closed and a constitutive model for the stress is necessary. 
We will present the model employed in this study in the next section. 

2.2. Constitutive model 

We briefly present the constitutive model in its general form for com- 
pleteness. Readers are referred to |T| for details on the derivation and char- 
acteristics of the model. In this model, the stress is expressed as 

<r = p\ - p V S + a 3 p(AS + SA - | (A : S)l) + a A p(fk - -(A : S)S), (3) 

where p is the pressure or mean normal stress, S := S/ yj ^D T : D with S and D 

as the deviatoric and full strain rate tensors, respectively, i.e. S = D — |tr(D)l, 
7] is a scalar macroscopic friction coefficient and I is the identity tensor. The 
last two terms in Eq. [3] model the normal stress differences (NSD) in simple 
shear flow, without which the model predicts equal normal stresses for such 
flow. The pressure and the macroscopic friction coefficient are linked to 
coordination number Z 2 and fabric tensor A using 

pd/k = ( ai + a 2 \A\)(Z 2 -Z c ) a , (4) 

and 

rj = h + b 2 A : S, (5) 

respectively, where d is the particle diameter, k is the particle normal stiff- 
ness, and Z c a critical coordination number characterizing jamming. In 
Eqs. [3]-[5j the coefficients, a,, bj (i = 1,2,3,4 and j = 1,2) and a, are 
material constants with ax j2 and a independent of particle friction coefficient 
yti and the rest as a function of \l. Their numerical values are given in Table [T] 
The fabric evolution is modeled as 

A = ciS + c 2 |D|A + c 3 (A:S)A, (6) 

where A = A + A-W-W-A, with W as the spin tensor, W = ^(Vv - (Vv) T ), 
and A denotes its material time derivative. The coordination-number evolu- 
tion is given as 

Z 2 = d 1 {A : S - X |S|) + d 2 \D\{f(<f>) - Z 2 ) + d 3 tr(D), (7) 
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Table 1: Numerical values of a, ax, a2 and X1-X3 in the fitting expression xie 22 ^ + X3, for 
the material parameters in the pressure and r\ equations, and the Zi evolution equation [T]. 
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-0.52 -2.8 


100 
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7.5 


0.5 



Table 2: Numerical values of the material parameters in the microstructure evolution 
equations pp. 



where \ = — (02 + \/c\ — 8ciCs)/2c 3 equals A : S for steady simple shear 
according to Eq. |6j f{4>) = Z c + fi\{4> — (fic) 132 and tr(D) is the trace of D. In 
Eqs.|6]and[7j the coefficients, Cj, di (i = 1, 2, 3) and f3j (j = 1, 2), are material 
constants modeled independent of /1, whose values are listed in Table [2] 



Using the assumptions for radial hopper flow in section 2.1 the strain-rate 



tensor is diagonal, with [D rr , D ee , D w \ = [ff, v f, v f] giving S = 
This together with the fabric Eq. [6] determines that the fabric tensor is also 
diagonal, i.e. A = [A rr , A ee , A w ], with A ee = A w . This means —\A rT = 
Ag g = A w since the fabric tensor is defined as being traceless. The fabric 
evolution equation is thus simplified to 



dA r y 
dr 



- 1 (2ci + V3c 2 A rr + 3c 3 A 2 rr ). (8) 



The boundary condition of A rr = at the top surface of radius r\ (see Fig. [T] 
(b) for its location) is applied for this equation, which is a reasonable assump- 
tion of the isotropic state due to random deposition and lack of substantial 
shear deformation at the top. 

Substituting the fabric and strain-rate tensors in Eq. [5j we obtain 

77 = 61 + y/3b 2 A rr . (9) 
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Similarly with the stress Eq. |3j we have the stress components 



a rr = p ( 1 - ^n^j , (10) 
1 



oee = o w = V ( 1 + -/if, 

where ^ is a bulk friction coefficient defined as the ratio of (Von Mises) shear 
stress to pressure and 

fi b = V3(r) - a 3 A rr ). (11) 

The bulk friction coefficient can thus be calculated from fabric and particle 
friction and need not be taken as a constant as in many plasticity theories. 
Defining the relation between the stress components as Ka rr = a ee = a w , 
we then get the stress ratio 

K=\±^. (12) 

It can be seen that the relation between the stress components is derived 
from the stress equation given the fabric and flow conditions. We do not 
invoke the Haar-von Karman hypothesis [3] to reach this result. 

With this relation for stress components, the momentum conservation 
given in Eq. [2] can be recast as 

d&rr = 2<PV 2 | 2(K-l)a Tr 
df f 5 f 

For this equation, we apply two boundary conditions, assuming a free-fall 
arch of radius r and stress- free top surface of radius r\. We note that the 
free-fall arch concept is echoed by the force-chain idea in recent granular 
material research. Such arch can be seen from the visualization of force 



chains in the DEM simulation of a flat-bottom silo [13]. Equation 13 



is 



closed by the constitutive relations Eqs. [8- 12 and specific solutions can be 



obtained applying the boundary conditions. The pressure and Z2 evolution 
equations are redundant for incompressible flow, as evidenced from the above 
derivation process. 
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3. Hopper flow predictions 



Solving the radial model obtained in section 2J2 for hopper flow of par- 
ticles with friction /i = 0.3, the stress and fabric radial distributions are 
obtained and shown in Fig. [2] The present theory predicts that the normal- 
stress ratio K increases along the flow direction (see Fig. [2] (a)), which nat- 
urally arises from the variation of the fabric (i.e. A rr , shown in Fig. [2] (b)). 
The value A rr increases along the flow, indicating microstructure anisotropy 
increased due to large deformation. As prescribed in Eq. |9j the increase 
of A rr results in the increase in rj along the flow direction, which in turn 
increases the K values according to Eq. [12} The NSD effect further cou- 



ples the A rr variation to K, reflected in Eq. 11 As seen in Fig. 2 (a), the 



formulation with the NSD leads to smaller K values, 1.8 at the bottom, com- 
paring to K = 2.1 calculated without the NSD. In contrast, a constant K is 
prescribed by the Mohr-Coulomb yield criterion together with the incipient- 
yield-everywhere assumption, which was employed in many previous hopper- 
flow analyses [2U El more references therein], and simulations [26]. The 
essential difference in the present constitutive model is that the K value, a 
bulk level property, need not be specified a priori; but is calculated from the 
particle-level property, fi, and the microstructure evolution. Another dis- 
tinction has to be made from the constant K used in pressure analysis for a 
static granular column, which appears to be a reasonable assumption in the 
static condition, as seen in Janssen's analysis [3] and verified more recently 
by an elasticity theory [27]. This comparison highlights the fact that the 
stress state in the flow is different from that in the static condition. 

To compare these continuum results with DEM ones, we calculate an 
internal friction coefficient /, using the same definition as in a DEM study 
of two-dimensional hopper flow [12] . This coefficient is related to the stress 
ratio K as / = The value of / increases from 0.121 at the top to 

0.382 (0.300) at the bottom for the case without (with) the NSD as plotted 
in Fig[2^c). This variation agrees with that of 0.150-0.389 from the DEM 
simulation for a hopper with a = 30° discharging polydisperse particles of 
/i = 0.3 [12]. We note that it is reasonable to compare to the polydisperse 
system simulated as the monodisperse system formed hexagonal crystalline 
structures [12] . which is not consistant with the physical basis of the current 
constitutive model and not realistic for practical hopper flows. This trend 
of internal friction increasing from top to bottom is also supported by DEM 
simulations of cylindrical flat-bottom silos [TT) IT3"] , where the ratios between 
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the reported shear and normal stresses varied in a similar fashion. How- 
ever, no quantitative comparison can be made due to the lack of principal 
stress data and the different geometry used in [TU [13]. With support from 
the DEM simulations, our results confirm that the new constitutive model 
has advantages in predicting realistic characteristics of stress in the conical 
hopper flow. 

Predictions of normal stress a„ from the present formulations with and 
without the NSD effect and from the HGT are compared in Fig. [2] (d). The 
HGT uses a constant K equal to the mean of the K values from the calcu- 
lation with NSD. The profiles of a rr along the radial distance f are similar, 
with the HGT result being the largest and reaching the peak value at the 
lowest position. The differences in the normal stress profiles also lead to 
different predictions of the flow rate. The values of V from the calculations 
without NSD, with NSD and using HGT are 1.550, 1.867 and 1.992, respec- 
tively, which shows that larger a rr results in a higher flow rate. To put this 
flow rate in perspective to commonly measured discharge rates, we show that 
the Beverloo equation in the form of W = Cp<J)>Jg(D — kd) 5 [281 E]> where 
D is the aperture diameter and k is of order unity, can be related to V as 
sin a). The C measured in experiments takes a value close 
to 0.58 [4j. However, according to the predictions of V, C is at least 1.22, 
about twice of the common value for a a = 30° hopper. This over-prediction 
should mostly be ascribed to the radial analysis itself, not the constitutive 
model though. One major source of error is the assumption of smooth wall 
in this analysis. Because of this gross over-prediction, we shall not assess the 
capability of the three different analyses regarding discharge rate prediction. 
As a summary of the NSD effect, with NSD the predictions of the normal 
stress a rr and the flow rate are higher, but the stress ratio K is lower, than 
without NSD. 

Another important functionality of the constitutive model used in this 
study is to directly probe the particle-friction effect on flow behavior. We 
have performed the same hopper-flow analysis for granular materials with 
particle friction /x = 0.1 and /x = 1.0. The results for K and plotted 
in Figs. [3] (a)-(b), and (c)-(d), respectively, of which (a)-(c) are for /x = 0.1 
and (b)-(d) for /x = 1.0. The fabric is modeled as independent of /x and 
thus not repeated in this figure. Comparing Figs. [3] (a) and (b) together 
with Fig. [2] (a), it can be seen that increasing /x results in higher K values. 
The difference between predictions with and without NSD is also enlarged. 
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No NSD 

With NSD 




/ 



(c) (d) 

Figure 2: Radial variation of (a)stress ratio, (b) component of the fabric tensor, (c) macro- 
scopic friction coefficient and (d) normal stress, plotted against the dimensionless distance 
in the vertical axis. The material has particle friction = 0.3. 
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Similar trend is observed for the internal friction coefficient / (not shown), 
which is consistent with engineering observations although no quantitative 
comparison can be made at the moment. Studying the a rr in Figs. [3] (c)-(d) 
and Fig. [2] (d), it can be concluded that higher [i leads to lower normal stress 
and does not alter the effect of NSD. From the flow-rate analysis above, it can 
also be deduced that higher \i results in lower flow rates. Specifically, the V 
values from the calculations without NSD, with NSD, and using HGT equal 
2.161, 2.403 and 2.545, respectively, when \i = 0.1; and 1.365, 1.718 and 
1.83, respectively, when fi = 1.0. This effect of particle friction qualitatively 
matches practical experience and is supported by the DEM simulations of a 
cylindrical silo [16J, but again quantitative comparison cannot be made due 
to the lack of directly comparable data. The increase in friction also amplifies 
the difference between the flow rate predictions for cases with and without 
NSD from about 11% to 26%. However, the potential of the constitutive 
model in facilitating understanding of the particle-property influence on the 
stress distribution and flow rate can clearly be seen from the above results. 

4. Conclusions 

A radial-flow solution of the conical hopper flow has been carried out 
utilizing a recently-developed constitutive model [Tj. Results on stress, mi- 
crostructure and flow rate have been obtained by simultaneously solving the 
continuity, momentum and fabric evolution equations. It has been demon- 
strated that a varying stress ratio is predicted, as a result of the fabric evolu- 
tion. The particle friction has been shown to lead to higher stress ratios, but 
lower normal stress and flow rates as it increases. These results are supported 
qualitatively and quantitatively by DEM simulations and engineering obser- 
vations. The effects of considering the NSD, have been explored, showing 
that considering NSD results in higher normal stress and flow rate, but lower 
stress ratio. The difference can be as large as about 30%, indicating that the 
NSD in constitutive models warrants more attention in order to make more 
accurate predictions. 

The significance of the constitutive model has been demonstrated to lie in 
connecting particle-level properties to continuum-level computation through 
microstructure evolution, obviating the need for specifying bulk properties 
and enabling prediction of particle-property effects on flow. Admittedly, the 
radial analysis performed is too simplified to validate against other spatial 
measurements, or to be used as a predictive tool for hopper design. Fu- 
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Figure 3: Radial variation of stress ratio in (a) and (b), of normal stress in (c) and 
(d), plotted against the dimensionless distance in the vertical axis. The particle-friction 
coefficient is fi — 0.1 for (a) and (c), and fi = 1.0 for (b) and (d). 
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ture work has been planned on applying the constitutive model to multi- 
dimensional analysis of compressible flow and validating against salient hop- 
per flow behavior, such as shear localization. 
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